%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Code to reconstruct aggregate UR from counterfactual flows
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

close all
clear all

global jei jeu jui jue jie jiu

addpath ('ImpliedHazardRates')
addpath ('lib')
addpath ('lib\4Figures')

% Load hzd rates adjusted for demographics by DFM 
jei=xlsread(['Original.xlsx'],'JEI');
jeu=xlsread(['Original.xlsx'],'JEU');
jue=xlsread(['Original.xlsx'],'JUE');
jui=xlsread(['Original.xlsx'],'JUI');
jie=xlsread(['Original.xlsx'],'JIE');
jiu=xlsread(['Original.xlsx'],'JIU');


%LF shares
a=['w1624'; 'w2534'; 'w3544'; 'w4554'; 'w5564'; 'm1624'; 'm2534'; 'm3544'; 'm4554'; 'm5564'; '_6585'];

for index=1:length(a)
    data=xlsread(['FlowData\Flows\Flows_',num2str(a(index,:)),'.xlsx']);
    
    data=data(2:end-2,:);
    
    E(:,index)=data(:,7);
    U(:,index)=data(:,8);
    I(:,index)=data(:,9);
    
end

%UR BLS
u_bls=U./(U+E);

%LF shares
w=(E+U)./kron(sum(E+U,2),ones(1,size(E,2)));
wbar=mean(w);
%Pop shares
wpop=(E+U+I)./kron(sum(E+U+I,2),ones(1,size(E,2)));
wpopbar=mean(wpop);

% Published UR and LFPR
Ur_BLS=sum(U,2)./(sum(U+E,2));
Lr_BLS=sum(U+E,2)./(sum(U+E+I,2));


% Construct actual and counterfactual U rates:
jei0=jei;jui0=jui;jue0=jue;jeu0=jeu;jie0=jie;jiu0=jiu;
for kk=1:2
    jei=jei0;jui=jui0;jue=jue0;jeu=jeu0;jie=jie0;jiu=jiu0;
    
    if kk==2
        
        % Load counter-factual hzd rates 
        
                        jei=xlsread(['Original_adjusted.xlsx'],'JEI');
                        jeu=xlsread(['Original_adjusted.xlsx'],'JEU');
                        jue=xlsread(['Original_adjusted.xlsx'],'JUE');
                        jui=xlsread(['Original_adjusted.xlsx'],'JUI');
                        jie=xlsread(['Original_adjusted.xlsx'],'JIE');
                        jiu=xlsread(['Original_adjusted.xlsx'],'JIU');
         
    end
    
    %store factual or counter-factual flows:
    jui_sim(:,:,kk)=jui;
    jei_sim(:,:,kk)=jei;
    jie_sim(:,:,kk)=jie;
    jiu_sim(:,:,kk)=jiu;
    jue_sim(:,:,kk)=jue;
    jeu_sim(:,:,kk)=jeu;
    
    s=jei.*jiu+jie.*jeu+jiu.*jeu;
    f=jui.*jie+jiu.*jue+jie.*jue;
    o=jeu.*jui+jue.*jei+jui.*jei;
    
    Ur_ss{kk}=s./(s+f);
    Lr_ss{kk}=(s+f)./(s+f+o);
    
    Pop=U+E+I;
    for z=1:size(jei,2)
        %Now get the actual Ur and LFPr from SS values
        ul=lom_dis( [1-exp(-jei(:,z)) 1-exp(-jeu(:,z)) 1-exp(-jui(:,z)) 1-exp(-jue(:,z)) 1-exp(-jie(:,z)) 1-exp(-jiu(:,z))],Pop(:,z),U(:,z),E(:,z),I(:,z)); %code that iterates to get actual stocks
        Ur_temp(:,z)=ul(:,1);
        Lr_temp(:,z)=ul(:,2);
    end
    
    Ur{kk}=Ur_temp;
    Lr{kk}=Lr_temp;
    
    if kk==1
        % 1) Construct ACTUAL (flow-based) U and LFP rates:
        % Construct Agg LFP rate
        Lr_agg{kk}=sum(Lr{kk}.*wpop,2);
        % Construct Agg U rate
        Ur_agg{kk}=1./Lr_agg{kk}.*sum(Ur{kk}.*Lr{kk}.*wpop,2);
    end
    if kk==2
        % 2) Construct COUNTER-FACTUAL (flow-based) U and LFP rates:
        % Construct Agg LFP rate (use wpopbar, because keep demog constant)
        Lr_agg{kk}=sum(Lr{kk}.*kron(wpopbar,ones(length(Lr{kk}),1)),2);
        % Construct Agg U rate (use wpopbar, because keep demog constant)
        Ur_agg{kk}=1./Lr_agg{kk}.*sum(Ur{kk}.*Lr{kk}.*kron(wpopbar,ones(length(Lr{kk}),1)),2);
    end
end

% 
% LFSS-based demographic-adjused UR
Ur_agg_LFSS=sum(Ur{1}.*kron(wbar,ones(length(Lr{1}),1)),2);
% 
% LFSS-based demographic-driven UR
Ubar=mean(U./(U+E),1);
epsw_lf=kron(Ubar,ones(length(w),1)).*(w-kron(wbar,ones(length(jue),1)));
demog_lf=sum(epsw_lf,2);

%Pop-based demogr-driven UR
Urbar_agg=mean(Ur_BLS);
LFPRbar_agg=mean(Lr_BLS);
LFPRbar=mean((U+E)./(U+E+I),1);
lil=kron(LFPRbar./LFPRbar_agg,ones(length(jue),1));
ui_u=kron(Ubar-kron(Urbar_agg,ones(1,size(jei,2))),ones(length(jue),1));
epsw_pop=ui_u.*lil.*(wpop-kron(wpopbar,ones(length(jue),1)));
demog_pop=sum(epsw_pop,2);

%% FIGURES
% Plot figure 4 and 5 (actual and cftl UR)
Cftl_UrAgg

%Plot figure 8 (some counterfactual flows with actual flows) 
Cftl_flows

% Illustrate implications of counter-factuals for different subgroups
% PRIME-AGE WOMEN:
Cftl_women
% YOUNG:
Cftl_young  

%data for figures 6 and 7
qy(mq([trend_w trend_yg]));


